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We  present  a  simulation  tool  to  study  fluid  mixtures  that  are  simultaneously  chemically  reacting  and 
adsorbing  in  a  porous  material.  The  method  is  a  combination  of  the  reaction  ensemble  Monte  Carlo 
method  and  the  dual  control  volume  grand  canonical  molecular  dynamics  technique.  The  method, 
termed  the  dual  control  cell  reaction  ensemble  molecular  dynamics  method,  allows  for  the 
calculation  of  both  equilibrium  and  nonequilibrium  transport  properties  in  porous  materials  such  as 
diffusion  coefficients,  permeability,  and  mass  flux.  Control  cells,  which  are  in  direct  physical  contact 
with  the  porous  solid,  are  used  to  maintain  the  desired  reaction  and  flow  conditions  for  the  system. 

The  simulation  setup  closely  mimics  an  actual  experimental  system  in  which  the  thermodynamic 
and  flow  parameters  are  precisely  controlled.  We  present  an  application  of  the  method  to  the  dry 
reforming  of  methane  reaction  within  a  nanoscale  reactor  model  in  the  presence  of  a  semipermeable 
membrane  that  was  modeled  as  a  porous  material  similar  to  silicalite.  We  studied  the  effects  of  the 
membrane  structure  and  porosity  on  the  reaction  species  permeability  by  considering  three  different 
membrane  models.  We  also  studied  the  effects  of  an  imposed  pressure  gradient  across  the  membrane 
on  the  mass  flux  of  the  reaction  species.  Conversion  of  syngas  (H2/CO)  increased  significantly  in 
all  the  nanoscale  membrane  reactor  models  considered.  A  brief  discussion  of  further  potential 
applications  is  also  presented.  ©  2004  American  Institute  of  Physics.  [DOI:  10.1063/1.1782031] 


I.  INTRODUCTION 

With  the  rapid  growth  of  nanotechnology  and  the  inven¬ 
tion  of  various  nanomaterials,  some  of  these  materials  have 
been  proposed  as  vehicles  for  nanochemical  devices  such  as 
nanoscale  reactors  and  nanoscale  membrane  reactors.1  How¬ 
ever,  development  of  these  applications  is  impossible  with¬ 
out  fundamental  knowledge  of  reaction,  adsorption,  and 
transport  mechanisms  in  the  nanoporous  materials. 

It  is  well  established  that  confinement  brings  about  dras¬ 
tic  changes  in  the  thermodynamics  properties  of  fluids  such 
as  narrowing  of  the  coexistence  curve,  lowering  of  the  pore 
critical  temperature,  or  increasing  the  average  pore  densities; 
for  a  comprehensive  review  see  Ref.  2.  Confinement  also 
influences  chemical  reaction  equilibrium.  For  example,  gen¬ 
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erally  the  pore  phase  has  a  higher  density  than  the  corre¬ 
sponding  bulk  phase;  this  results  in  an  increase  in  yield  for 
reactions  in  which  there  is  a  decrease  in  the  total  number  of 
moles  (Le  Chatellier’s  principle).  Further,  some  components 
of  the  reaction  mixture  are  selectively  adsorbed  on  the  solid 
surfaces,  subsequently  affecting  the  reaction  equilibrium.  Fi¬ 
nally,  molecular  orientations  can  be  strongly  influenced  by 
proximity  to  a  solid  surface  which  also  can  shift  the  reaction 
equilibrium  from  the  bulk  phase  equilibrium.  The  effects  of 
confinement  on  reaction  equilibria  were  already  studied  by 
Turner  et  al?  for  several  realistic,  reversible  reactions  in 
carbon  micropores  and  carbon  nanotubes,  and  by  Borowko 
et  al.6,7  for  model,  reversible  reactions  in  slitlike  pores. 

These  authors  employed  the  reaction  ensemble  Monte 
Carlo  (REMC)  simulation  technique8-10  that  enables  us  to 
directly  simulate  equilibrium  properties  of  chemically  react¬ 
ing  systems.  The  method  requires  only  a  knowledge  of  the 
reaction  species  intermolecular  potentials  and  their  ideal-gas 
properties,  in  addition  to  specification  of  the  system  stoichi- 

©  2004  American  Institute  of  Physics 


4901 


Downloaded  16  Aug  2006  to  128.63.163.77.  Redistribution  subject  to  AIP  license  or  copyright,  see  http://jcp.aip.org/jcp/copyright.jsp 


4902  J.  Chem.  Phys.,  Vol.  121,  No.  10,  8  September  2004 

ometry  and  thermodynamic  constraints.  Recently,  we  pro¬ 
posed  a  combination  of  the  REMC  method  with  the  molecu¬ 
lar  dynamics  (MD)  technique.11  The  method,  termed  the 
reaction  ensemble  molecular  dynamics  (RxMD)  method, 
uses  a  combination  of  stochastic  and  dynamic  simulation 
steps,  allowing  for  the  simulation  of  both  thermodynamic 
and  transport  properties.  The  method  couples  a  MD  system 
(dynamic  cell)  to  a  reaction  mixture  reservoir  (control  cell) 
that  is  formulated  upon  the  REMC  method.  Thermodynamic 
and  transport  properties  are  calculated  in  the  dynamic  cell  by 
using  a  canonical  MD  simulation  method.  The  RxMD 
method  is  analogous  to  the  grand  canonical  molecular  dy¬ 
namics  (GCMD)  technique  of  Papadopoulou  et  al. 12  The  ac¬ 
curacy  and  stability  of  the  RxMD  method  was  assessed  by 
considering  the  ammonia  synthesis  reaction,  N2  +  3H2 
^2NH3. 

Confinement  also  influences  the  transport  properties  of 
fluid  particles  inside  the  nanoporous  materials.  Physical 
space  restrictions  based  on  the  fluid  particle  size  or  geometry 
may  limit  flow  of  particles  through  particular  pores  in  the 
material.  Furthermore,  attraction  of  the  pore  surface,  i.e., 
physisorption  may  play  a  critical  role.  Even  further  compli¬ 
cating  the  matters,  fluid  particles  may  chemisorb.  In  addition 
to  phenomena  occurring  between  fluid  particles  and  the  pore 
surface,  behavior  in  the  porous  material  becomes  increas¬ 
ingly  more  complex  if  chemical  reactions  occur  between 
fluid  particles. 

Confinement  contributes  significantly  to  the  thermody¬ 
namic  and  transport  properties  if  reactions  and  separation 
occur  simultaneously  in  the  nanoporous  materials.  Consider, 
for  example,  two  voids  in  a  nanoporous  material  separated 
by  a  semipermeable  membrane  with  the  elementary  revers¬ 
ible  reaction 

A+B^C  +  D  (1) 

occurring  in  one  void  only.  If  the  membrane  is  permeable  to 
product  D  only  and  the  partial  pressure  of  D  on  the  permeate 
side  of  the  membrane  is  lower  than  the  partial  pressure  of  D 
on  the  reaction  side  of  the  membrane  then  the  separation 
process  produces  a  flow  of  product  D  across  the  membrane. 
The  removal  of  D  from  the  “reaction”  void  via  the  separa¬ 
tion  function  of  the  membrane  has  the  following  effects  on 
the  reaction  shown  in  Eq.  (1):  (i)  the  reaction  equilibrium 
condition  is  shifted  to  the  right  resulting  in  a  higher  conver¬ 
sion  of  reactants  A  and  B  to  products  C  and  D ;  (ii)  if  there 
is  an  undesirable  side  reaction,  e.g.,  that  fouls  a  catalyst  by  a 
side  product,  such  as 

B  +  D^±E  (2) 

occurring  in  the  reaction  void,  then  the  separation  of  product 
D  from  the  reaction  mixture  reduces  the  amount  of  reactant 
B  to  the  side  reaction,  increasing  the  selectivity  of  conver¬ 
sion  to  product  C  (or  Z)).13 

In  this  work,  we  propose  a  nonequilibrium  MD  method 
for  the  simulation  of  combined  reaction  and  adsorption 
mechanisms  in  porous  materials.  This  method  is  a  combina¬ 
tion  of  the  REMC  and  RxMD  methods,  and  the  dual  control 
volume  GCMD  (DCV-GCMD)  technique.  The  REMC  and 
RxMD  methods  predict  the  physical  effects  on  chemical  re- 


Li’sal  et  at. 

action  equilibria  from  such  influences  as  solvation  and  con¬ 
finement.  In  addition,  the  RxMD  method  predicts  the  effects 
of  nonideal  environments  on  the  molecular  transport  proper¬ 
ties  of  chemically  reacting  mixtures.  The  DCV-GCMD 
method  was  independently  proposed  by  Heffelfinger  and 
van  Swol,14  and  MacElroy,15  and  was  extensively  employed 
for  studies  of  transport  phenomena  in  confined  systems;  see, 
e.g..  Refs.  16  and  17,  and  references  therein. 

We  term  the  nonequilibrium  simulation  method  the  dual 
control  cell  reaction  ensemble  molecular  dynamics  (DCC- 
RxMD)  method.  The  DCC-RxMD  method  allows  for  the 
simulation  of  the  thermodynamic  and  transport  properties  of 
fluids  in  porous  materials  where  reaction  and  adsorption  is 
occurring  simultaneously.  Fluid  particles  move  through  the 
simulation  cell  via  MD,18  chemical  reactions  occur  via  the 
REMC  method,  and  fluid  transport  (occurring  as  a  result  of 
an  imposed  pressure  gradient  between  adjacent  pores)  via  the 
grand  canonical  Monte  Carlo  (GCMC)  method.19 

We  apply  the  DCC-RxMD  method  to  study  reactions 
and  separations  in  a  model  of  a  nanoscale  membrane  reactor 
for  dry  reforming  of  methane.  Dry  reforming  of  methane  is 
an  important  industrial  process  for  producing  syngas 
(H2/CO)  that  utilizes  membrane  technology.20  We  examine 
the  effects  of  the  membrane  structure  and  porosity  on  the 
reaction  species  permeability  by  considering  three  different 
membrane  models.  We  also  examine  the  effects  of  an  im¬ 
posed  pressure  gradient  across  the  membrane  on  the  mass 
flux  of  the  reaction  species  and  reaction  conversion. 

The  paper  is  organized  as  follows:  Derivation  and  gen¬ 
eral  computational  considerations  of  the  DCC-RxMD  meth¬ 
odology  are  presented  in  Sec.  II.  Section  III  describes  an 
illustrative  application  of  the  DCC-RxMD  method  for  simu¬ 
lation  of  a  multicomponent  system  in  which  both  a  chemical 
reaction  (dry  reforming  of  methane)  and  physisorption  is  oc¬ 
curring  within  a  nanoporous  solid.  In  addition,  the  technical 
details  of  this  application  of  the  DCC-RxMD  methodology 
are  given.  Section  IV  presents  and  discusses  DCC-RxMD 
results  together  with  DCV-GCMD  simulations  for  both  the 
pure  fluids  and  equimolar  mixtures  of  the  reaction  compo¬ 
nents  (CH4,  C02,  Hi,  CO).  Section  V  contains  our  conclu¬ 
sions  including  a  brief  discussion  of  further  potential  appli¬ 
cations  of  the  DCC-RxMD  method. 

II.  SIMULATION  METHODOLOGY 

A  molecular  simulation  method  to  study  reactions  and 
adsorption  in  porous  materials  is  presented.  The  method  uses 
a  combination  of  stochastic  and  dynamic  simulation  steps, 
allowing  for  the  simulation  of  both  thermodynamic  and 
transport  properties.  The  method  couples  a  MD  system  (dy¬ 
namic  cell)  to  reaction  and  nonreaction  mixture  reservoirs 
(control  cells)  that  are  formulated  upon  the  REMC 
method8-10  and  GCMC  method,19  hence  the  term  dual  con¬ 
trol  cell  reaction  ensemble  molecular  dynamics  (DCC- 
RxMD)  method.  The  control  cells  are  in  direct  contact  with 
the  dynamic  cell  and  the  particles  are  able  to  move  freely 
between  the  cells.  Transport  properties  are  calculated  in  the 
dynamic  cell  by  using  a  MD  simulation  method. 18  REMC 
forward  and  reverse  reaction  steps,  and  GCMC  particle  in¬ 
sertion  and  deletion  steps  are  performed  in  the  control  cells 
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FIG.  1.  Schematic  of  a  nanoscale 
membrane  reactor  model  made  up  of  a 
reaction  void  and  a  transport  void 
separated  by  a  semipermeable  mem¬ 
brane.  A  model  reaction  2 A^B  takes 
place  in  the  reaction  void  while  com¬ 
ponent  B  is  separated  via  the  semiper¬ 
meable  membrane.  Enlargement 
shows  schematic  of  the  corresponding 
DCC-RxMD  simulation  setup.  Peri¬ 
odic  boundary  conditions,  applied  in 
both  the  y  and  z  directions,  are  omit¬ 
ted  in  the  x  direction  due  to  the  pres¬ 
ence  of  repulsive  confining  walls  (in 
the  yz  plane)  at  each  end  of  the  simu¬ 
lation  setup. 


only,  while  MD  steps  are  performed  in  both  the  dynamic  cell 
and  the  control  cells.  The  control  cells,  which  act  as  sink  and 
source  reservoirs,  are  maintained  at  explicitly  defined  and 
controlled  gradient  conditions,  e.g.,  a  pressure  gradient.  The 
DCC-RxMD  method  is  analogous  to  the  DCV-GCMD 
method; 14,15  all  of  which  simulate  conditions  that  directly 
relate  to  real,  open  systems.  Computational  details  will  be 
addressed  below  but  first  we  illustrate  the  overall  setup  of  the 
method. 

Consider  a  nanoporous  material  containing  a  system  of 
“nanomembrane  reactors.”  Each  nanomembrane  reactor  is 
made  up  of  two  adjacent  voids  separated  by  a  semipermeable 
membrane.  One  of  the  voids  acts  as  a  “nanoreactor”  in 
which,  for  example,  the  model  reaction 

2  A^B  (3) 

is  occurring.  We  term  this  void  the  reaction  void.  The  ac¬ 
companying  void  in  the  nanomembrane  reactor,  termed  the 
transport  void,  maintains  a  pressure  gradient  across  the 
membrane.  We  further  assume  that  the  membrane  separating 
the  reaction  and  transport  voids  is  not  permeable  to  all  mo¬ 
lecular  components.  For  example,  in  the  model  system  of  Eq. 
(3),  suppose  the  membrane  is  permeable  to  component  B 
only.  For  such  a  case,  a  difference  in  the  partial  pressures  of 
component  B  in  the  reaction  and  transport  voids  results  in  a 
flux  of  component  B  through  the  membrane.  A  schematic  of 
the  nanomembrane  reactor  model  made  up  of  the  reaction 
and  transport  voids  separated  by  a  semipermeable  membrane 
is  presented  in  Fig.  1.  The  nanomembrane  reactor  model  is 
analogous  to  the  slitpore  model21  which  assumes  that  the 
nanomembrane  reactor  model  is  taken  as  representative  of  all 
nanomembrane  reactors  in  the  porous  material.  The  effect  of 
pore  connectivity  can  be  readily  involved  by  considering  a 
combination  of  connected  nanomembrane  reactor  models. 

To  mimic  reaction  and  adsorption  mechanisms  in  such  a 
nanoscale  membrane  reactor  model,  we  consider  a  DCC- 
RxMD  simulation  box  as  shown  in  the  enlargement  of  Fig.  1. 


A  membrane  of  thickness  S—2Lix  is  placed  at  the  center  of 
the  simulation  box.  The  left-hand  side  of  the  simulation  box 
(from  ~L\X  to  —  L4x)  with  control  cell  (CC)  A  (from  —  L2x 
to  —Lyx)  corresponds  to  a  reaction  void  while  the  right-hand 
side  of  the  simulation  box  (from  L]x  to  L4x)  with  CC  B 
(from  L2x  to  L2x)  corresponds  to  a  transport  void.  Box 
lengths  in  the  y  and  z  directions  are  denoted  by  Ly  and  Lz , 
respectively.  For  simplicity,  we  considered  that  both  CC  A 
and  CC  B  are  of  the  same  size  and  are  positioned  symmetri¬ 
cally  with  respect  to  the  yz  plane  at  x  =  0.  We  impose  peri¬ 
odic  boundary  conditions  in  both  the  y  and  z  directions  since 
we  assume  that  the  sizes  of  the  reaction  and  transport  voids 
in  these  directions  are  much  larger  than  the  membrane  thick¬ 
ness.  To  limit  the  size  of  the  simulation  setup  in  the  x  direc¬ 
tion,  repulsive  confining  walls  are  placed  at  the  end  of  the 
simulation  box.  The  size  of  the  simulation  box  in  the  x  di¬ 
rection  must  be  sufficiently  large  to  ensure  that  these 
confining-wall  boundaries  have  negligible  effect  on  the  cal¬ 
culated  fluid  properties.  The  alternative  to  including 
confining-wall  boundaries  is  to  allow  the  fluid  particles  to 
leave  the  simulation  box  at  the  boundaries.  Use  of  either  the 
confining-wall  boundaries  or  open-end  boundaries  give  the 
same  simulation  results.15'22  An  alternative  scenario  is  pos¬ 
sible  here.  For  example,  one  can  consider  a  DCC-RxMD 
simulation  box  with  two  confining  walls  located  in  the  z 
direction  at  a  distance  L,  apart  without  confining  walls  in  the 
x  direction.22  In  this  scenario,  confinement  effects  in  the  z 
direction  are  also  considered  important. 

In  either  scenario,  the  DCC-RxMD  simulation  proceeds 
as  follows.  After  n  MD  MD  steps,  the  system  is  frozen,  i.e., 
particle  positions  are  held  fixed,  and  we  perform  «remc  for- 
ward  and  reverse  reaction  steps  in  CC  A,  and  «gcmc  particle 
creation  and  destruction  steps  in  CC  B .  Both  the  REMC  and 
GCMC  algorithms  require  insertion  of  particles  into  the  CCs. 
The  velocities  for  such  particles  are  assigned  from  a 
Maxwell-Boltzmann  distribution  corresponding  to  the  speci- 
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fled  system  temperature.11'12  Values  of  nMD,  Hremc,  and 
n  gcmc  must  be  chosen  appropriately  to  maintain  reaction 
equilibrium  in  CC  A,  constant  chemical  potentials  in  CC  B , 
and  reasonable  transport  rates  at  the  boundaries  between  the 
CCs  and  the  membrane  region  (the  dynamic  cell). 


A.  Isothermal  molecular  dynamics 


Trajectories  of  fluid  particles  within  the  entire  DCC- 
RxMD  simulation  volume  are  generated  by  the  MD  simula¬ 
tion  method.18  The  imposed  pressure  gradient  across  the 
membrane  (a  difference  between  pressures  in  CC  A  and  CC 
B )  produces  heat  that  must  be  removed  to  maintain  isother¬ 
mal  transport.  Removal  of  the  heat  can  be  achieved  by  keep¬ 
ing  the  kinetic  energy  of  the  system  Ek  constant  during  MD 
steps.23  Ek  is  defined  in  terms  of  Ap,  ;=p, ,  — p;,  the  relative 
momenta  of  particles  with  respect  to  the  streaming  momen¬ 
tum  per  particle  for  each  of  components  l.  Pi,  evaluated  as 


1  1 

P/=ttE  Pm 

W/  ;=i 


(4) 


where  Ar  is  the  time  step  and 


1 

=  77  2  f;,/- 
A/;=i 


(ID 


Ap  u  ( t )  is  then  used  to  compute  an  instantaneous  tempera¬ 
ture  of  the  system  T(  t )  from  the  relation 


jkB  T[t) 


1  1  ' 

~2  -2  (Apr;)2, 

2  i=\  mti=i 


(12) 


where  F=  3(E[=1./V/)  —  1  is  the  total  number  of  degrees  of 
freedom  for  the  system  and  kH  is  the  Boltzmann’s  constant. 
T{t)  and  the  specified  system  temperature  T  are  used  to 
evaluate  the  constant  ft  defined  as 

(13) 

where  ft  is  related  to  the  constant  a  [defined  by  Eqs.  (8)  and 
(9)]  via 


In  Eq.  (4),  /V)  is  the  number  of  particles  of  component  /  in 
the  system  and  p,  j  is  the  momentum  of  particle  i  for  com¬ 
ponent  /.  The  kinetic  energy  constraint  for  the  system  will  be 
satisfied24  if 


dEk 

dt 


1  d 

2  dt 


C 


m,i=  i 


(5) 


where  t  is  the  time,  m/  is  the  mass  of  component  /,  and  c  is 
the  total  number  of  components  in  the  system.  A  set  of  equa¬ 
tions  of  motion  for  the  MD  subject  to  the  kinetic  energy 
constraint  given  by  Eq.  (5)  is25 


dru  _  pM 

dt  nil  ’ 


(6) 


dPi,i 

dt 


=  f;,/-«(ApM), 


(7) 


where  r,  /  is  the  position  of  particle  i  for  component  /  and  f,  / 
is  the  force  on  such  a  particle,  a  is  computed  from  Eq.  (5) 


c  j  n, 

E  —  E  (f;,r  Ap,-  /) 

1  =  l  nii  i=  l 

E  —  2  ( Ap,  / -  Ap, ;) 

/=  i  ni/  i=  i 


1 

2 


C 


m-i  i=  i 


(f/./- Ap,,/) 


Ek 


(8) 


(9) 


The  equations  of  motion  for  the  isothermal  MD  method, 
Eqs.  (6)  and  (7),  were  solved  by  the  “leap-frog”  version  of 
the  Verlet  algorithm.26  The  algorithm  starts  with  Eq.  (7)  by 
evaluating  the  unconstrained  relative  momenta  Ap'"/1  at 
time  t. 


Ap“j(t)  =  Ap,/(?-At/2)+  y[fM(t)-f,(f)],  (10) 


13  = 


1 

1  +  a(At/2) ' 


(14) 


Having  determined  ft,  we  are  able  to  obtain  p,  /  at  time 
t  +  At/2  as 

Pi  i(t  +  At/2)  =  (2ft—  1  )p(j/(f  —  Af/2)  +2(/3—  1  )p/(f  —  Af/2) 

+At[ftfu(t)+(\-ft)f,m  as) 

Finally,  r, ,  at  time  t  +  At  is  obtained  from  Eq.  (6)  as 


rM(f  +  Af) 


r i,i(t)  +  At 


P;  j(t  +  At/2) 

mi 


(16) 


B.  Reaction  ensemble  Monte  Carlo 

The  REMC  method8-10  is  a  powerful  simulation  tool  for 
studying  chemically  reacting  mixtures.  The  method  only  re¬ 
quires  inputting  the  intermolecular  potentials  and  the  ideal- 
gas  properties  for  the  reaction  species  that  are  present.  Most 
notably,  the  method  does  not  require  a  reactive-type  potential 
that  mimics  bond  breakage  and  formation.27,28  The  REMC 
method  predicts  the  shift  in  equilibria  of  an  ideal-gas  phase 
reaction  due  to  nonideal  conditions  such  as  high  temperature 
and  high  pressure  as  well  as  nonideal  surrounding  environ¬ 
ments  such  as  solvents  and  pore  walls.  Reactions  are  simu¬ 
lated  by  performing  forward  and  reverse  reaction  steps  ac¬ 
cording  to  the  REMC  algorithm  which  guarantees  that  the 
reaction  equilibrium  criteria  for  a  set  of  R  linearly,  indepen¬ 
dent  chemical  reactions, 

C 

E  V  7/r,  =  0,  j  =  1,2, . .  .  ,R,  (17) 

/=i 

are  established.29  In  Eq.  (17),  n/7  is  the  stoichiometric  coef¬ 
ficient  of  component  /  in  chemical  reaction  j  and  /i;  is  its 
chemical  potential.  The  reaction  equilibrium  condition  for 
our  model  reaction,  Eq.  (3),  in  CC  A  (the  reaction  void)  is 
then 
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a4-2/4  =  0,  (18) 

where  fif  is  the  chemical  potential  of  component  l  in  CC  A . 
Forward  and  reverse  reaction  steps  are  accepted  with  prob¬ 
abilities 


min 


1, 


r  i\M-i) 


vA 


tf»+i 


-exp 


A  UA 

kRT 


and 


min 


1, 


V4 


r  (ivl+i)«+2) 


-exp 


A  UA 

knT 


(19) 


(20) 


respectively.  In  Eqs.  (19)  and  (20),  V A  is  the  volume  of  CC 
A,  T  is  the  ideal-gas  quantity  defined  as 


r= 


knT 


~K, 


(21) 


where  P°  is  the  standard  state  pressure  of  1  bar,  K  is  the 
equilibrium  constant,29  NA  is  the  number  of  particles  of  com¬ 
ponent  I  in  CC  A,  and  A U A  is  the  change  in  the  configura¬ 
tional  energy  U  due  to  forward  and  reverse  reaction  attempts 
in  CC  A . 


C.  Imposed  pressure  gradient  (GCMC) 


In  order  to  maintain  a  flux  of  particles  through  the  mem¬ 
brane,  we  impose  a  pressure  gradient  across  it.  The  pressure 
gradient  is  indirectly  controlled  by  performing  GCMC  par¬ 
ticle  insertion  and  deletion  steps  in  CC  B  (the  transport  void) 
only.  According  to  the  GCMC  algorithm,19  the  chemical  po¬ 
tential  of  component  /  in  CC  B ,  /if ,  is  chosen  and  the  cre¬ 
ation  and  destruction  of  a  particle  of  component  l  is  accepted 
with  probabilities 


min 


,  VBq,(T ) 


A  UBV 

U.rM 

kBT) 

and 


min 


1, 


N, 


VBq,(T) 


exp 


[  I 

"ST 

<1 

[  t,rlexp( 

kBT  1 

(22) 


(23) 


respectively.  In  Eqs.  (22)  and  (23),  VB  is  the  volume  of  CC 
B,  qj(T )  is  the  internal  contribution  of  the  partition  function 
for  component  l,  /\'f  is  the  number  of  particles  of  component 
l  in  CC  B,  and  A  UB  is  the  change  in  the  configurational 
energy  U  due  to  creation  and  destruction  attempts  in  CC  B. 


III.  ILLUSTRATIVE  APPLICATION:  A  NANOSCALE 
MEMBRANE  REACTOR  MODEL 


details  of  the  reaction  species  models,  membrane  models, 
and  further  computational  information  specific  to  such  the 
DCC-RxMD  application. 


A.  Reaction  species  potential  models 

We  consider  a  nanoscale  membrane  reactor  model  (de¬ 
scribed  in  the  preceding  section,  see  also  Fig.  1)  in  which  the 
dry  reforming  reaction 


CH4  +  C02^2H2  +  2C0 


(24) 


occurs  in  the  reaction  void  (CC  A)  and  the  membrane  is 
(preferentially)  permeable  to  H2 .  A  difference  in  the  H2  par¬ 
tial  pressures  between  the  reaction  and  transport  voids,  i.e., 
between  CC  A  and  CC  B  will  result  in  a  H2  flux  through  the 
membrane.  REMC  forward  and  reverse  reaction  steps  for  the 
reaction  shown  in  Eq.  (24)  are  accepted  with  probabilities 


hvlr 


NAchNAco2 

(<+!)(<  + 2  )(^Vco+  1  KNAC  0+2) 


(25) 


and 
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v\v  rCH+i)c+D 


-exp  I 


A  UA 

knT 


respectively,  where  F  is  defined  as 


dO  \  2 


r  = 


knT 


K 


(26) 


(27) 


and  K  is  evaluated  using  the  JANAF  Tables.10 

Components  of  the  reaction  system  CH4 ,  C02 ,  H2 ,  and 
CO  are  modeled  as  Lennard-Jones  (LJ)  spheres.  The  fluid- 
fluid  interactions  are  approximated  with  a  truncated- and  - 
shifted  LJ  (TS-LJ)  potential 


,TS-LJ/ 


(rii)  =  uL\rij)-u 


LJ 


(28) 


0’  0/  c .ab  , 
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\rc,abl 
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(31) 


We  applied  the  DCC-RxMD  method  to  simulate  a  mul¬ 
ticomponent  system  in  which  both  a  chemical  reaction  and 
physisorption  is  occurring  within  a  nanoporous  solid.  We 
consider  the  dry  reforming  of  methane  reaction  which  is  an 
important  industrial  process  for  producing  syngas  (H2/CO) 
from  CH4  and  C02 . 20  We  demonstrate  how  the  reaction  con¬ 
version  is  affected  by  separating  out  H2  from  the  reaction 
mixture  via  a  semipermeable  membrane.  Below  we  provide 


Here,  r;j  is  the  distance  between  particle  i  for  component  a 
and  particle  j  for  component  b,  e’s  and  cr’s  are  the  effective 
LJ  energy  and  size  parameters,  respectively,  and  rcab  is  the 
spherical  cutoff  distance  assigned  to  be  rc  ab  =  3.5<rab.  The 
s’s  and  cr’s  for  CH4,  C02,  and  H2  were  taken  from  Ref.  31, 
those  for  CO  were  taken  from  Ref.  32  and  are  listed  in  Table 
I.  The  cross-term  LJ  parameters  were  evaluated  using  the 
Lorentz-Berthelot  combining  rules,18 
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TABLE  I.  Effective  LJ  energy  (e)  and  size  (cr)  parameters  for  CH4 .  C02 , 
H2  (Ref.  31)  and  CO  (Ref.  32)  fluids,  and  membrane  particles  (Ref.  36). 


Component 

e  lkB 
(K) 

a 

(nm) 

ch4 

148.1 

0.3810 

co2 

225.3 

0.3794 

h2 

38.0 

0.290 

CO 

123.0 

0.3662 

Membrane 

82.0 

0.270 

£ab  ''l^aa^bb » 

(32) 

<Taa  +  u  bb 

(33) 

vab~  2  ■ 

Firouzi  et  al. 33  have  confirmed  that  these  LJ 

models  for 

CH4 ,  C02 ,  H2 ,  and  CO  produced  a  phase  diagram  and  ther¬ 
modynamic  properties  which  are  in  quantitative  agreement 
with  the  experimental  data,  under  both  subcritical  and  super¬ 
critical  conditions. 

Fluid  particles  of  component  /  interact  with  the  repulsive 
confining  walls  w  (see  Fig.  1)  via  the  TS-LJ  potential  with 
b,w,  (T,„,  and  rcJw  =  21/6cr/w.  For  simplicity,  we  use  slw 
=  e;/  and  (Thl.=  an  since  the  confining  walls  have  no  direct 
influence  on  the  transport  in  the  membrane  region  and  have 
only  small  effects  on  fluid  properties  in  the  portions  of  CCs 
adjacent  to  the  confining  walls. 

B.  Membrane  models 

In  order  to  increase  the  reaction  conversion  in  our  na¬ 
nomembrane  reactor  model  by  separating  out  a  particular 
product  from  the  mixture,  it  is  required  that  the  membrane 
separating  the  reaction  and  transport  voids  be  permeable  to 
that  particular  product  only.  For  the  dry  reforming  of  meth¬ 
ane  reaction  considered  here,  we  define  our  membrane  model 
to  be  permeable  to  H2  only  while  being  impermeable  to  the 
remaining  reaction  species,  namely,  CH4,  C02,  and  CO. 
Since  the  dry  reforming  reaction  is  typically  carried  out  at 
quite  high  temperatures,  in  general,  the  separation  of  mixture 
components  via  a  membrane  is  based  primarily  on  molecular 
sieving  caused  by  the  passage  of  smaller  molecules  of  the 
mixture  through  the  pores  while  the  larger  molecules  are 
obstructed.34  As  evident  in  Table  I,  H2  molecules  are  smaller 
than  the  CH4,  C02,  and  CO  molecules.  Hence,  a  character¬ 
istic  size  for  the  membrane  pores  should  be  greater  than  crH 
(to  allow  H2  to  permeate  through  the  membrane)  but  yet  less 
than  or  comparable  to  o-CH4tascrC02s“o'C0  (to  obstruct  CH4, 
C02,  and  CO  from  permeating  through  the  membrane). 

We  utilized  two  types  of  membrane  models.  The  first 
model,  proposed  by  Powles  el  al., 25,35  defines  membranes  by 
several  layers  of  LJ  particles  (characterized  by  the  energy 
and  size  parameters  em  and  crm ,  respectively)  with  a  dis¬ 
tance  between  layers  equal  to  2 1/6 crm  (which  corresponds  to 
the  LJ  potential  minimum).  The  particles  in  each  layer  are 
arranged  in  a  face-centered  cubic  (fee)  structure.  The  layers 
are  built  by  replicating  a  four-particle  fee  primitive  cell  in  the 
y  and  z  directions.  A  two-dimensional  schematic  of  a  layer  of 
particles  for  this  model,  termed  the  Powles’  membrane 


FIG.  2.  Schematic  of  a  fee  membrane  layer  for  the  Powles’  membrane 
model  discussed  in  text.  The  size  of  the  particles  making  up  the  membrane 
is  characterized  by  am  .  An  open  circle  of  diameter  indicates  the 

maximum  size  of  a  membrane  pore. 

model  (PM),  is  shown  in  Fig.  2.  Also  shown  is  the  maximum 
allowable  size  for  a  pore  in  the  membrane  structure. 

Chosen  values  for  em  and  crm  correspond  to  silicalite 
molecules36  and  are  listed  in  Table  I.  The  membrane  particles 
interact  with  fluid  particles  via  the  TS-LJ  potential,  Eqs.  (28) 
and  (29),  and  with  the  cross-term  LJ  parameters  evaluated 
from  the  Lorentz-Berthelot  combining  rules,  Eqs.  (32)  and 
(33).  Powles’  model  generates  membranes  with  well-defined 
structures  comprising  straight  channels  that  are  prototypical 
of  zeolitic  structures.  We  built  two  Powles’  membranes  and 
denote  them  PM1  and  PM2.  Both  models  consist  of  seven 
fee  layers  but  differ  by  their  values  of  ■  PM1  has 
=0.3253  nm  which  corresponds  to  a  membrane  number 
density  pm  =  21.728  nm~3  while  PM2  is  characterized  by 
=0.4104  nm  and  pm  =  16.635  nm” 3 .  Geometric  pore 
size  distributions37  (PSDs)  for  PM1  and  PM2  are  displayed 
in  Fig.  3.  As  expected.  Fig.  3  shows  PSDs  that  are  quite 
narrow  with  maximum  peaks  corresponding  to  d™*e  . 

The  second  type  of  membrane  model  used  in  this  study 
comprises  random  configurations  of  nonoverlapping  LJ 
spheres.  The  model,  termed  the  random  membrane  model 
(RM),  is  generated  from  a  canonical  Monte  Carlo  simulation 
of  hard  spheres  with  a  corresponding  diameter  of  crm  . 38  We 
built  RMs  with  various  PSDs  in  the  volume  2 LlxXLyXLz 
where  the  same  value  of  Llx  as  in  PM1  and  PM2  is  used.  We 
then  chose  one  RM  whose  PSD  peak  is  roughly  at  the  same 
position  as  the  PSD  peak  for  PM2.  The  chosen  RM  has  p„, 
=  15.156  nm~3  and  its  PSD  is  displayed  in  Fig.  3.  Figure  3 
shows  that  the  PSD  for  the  RM  is  substantially  broader  than 
that  of  PM2  due  to  the  random  arrangement  of  the  membrane 
particles. 

C.  Computational  details 

In  addition  to  the  DCC-RxMD  simulations,  we  also  car¬ 
ried  out  DCV-GCMD  simulations  for  pure  CH4,  C02,  H2, 
and  CO,  and  an  equimolar  mixture  of  CH4  /C02  /H2  /CO.  In 
the  DCV-GCMD  simulations  both  CCs  are  used  to  maintain 
constant  chemical  potentials  for  the  particular  components. 
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FIG.  3.  Computed  geometric  pore  size  distributions  (PSD)  for  the  three 

membrane  models:  PM1  ( - ),  PM2  ( - ),  and  RM  ( - ). 

Units  for  the  v  coordinates  are  arbitrary. 


Setting  different  values  for  the  chemical  potentials  in  CC  A 
and  CC  B  results  in  different  pressures  in  these  CCs.  The 
imposed  pressure  gradient  across  the  membrane  causes  a  flux 
of  fluid  particles  through  the  membrane.  Hence,  the  DCV- 
GCMD  simulations  enable  us  to  investigate  properties  of  the 
membranes  PM1,  PM2,  and  RM  independent  of  chemical 
reaction  behavior.  For  the  DCV-GCMD  and  DCC-RxMD 
simulation  boxes,  we  used  Llx=L2x=  0. 909066  nm  (Llt 
=  0.909  066  nm  corresponds  to  seven  membrane  layers  in 
PM1  and  PM2),  L2x=  18.288  nm,  L4x=  19.050  nm,  and  Ly 
=  LZ  =  9.525  nm.  Use  of  Lix=L2x  (see  Fig.  1)  minimizes 
possible  steric  hindrance  effects  at  the  entrance  of  the  mem¬ 
brane  pores.39 

After  some  preliminary  test  runs,  we  found  that  nMD 
=  10,  nREMC=300,  and  «gcmc=300  are  satisfactory  values 
to  properly  maintain  reaction  equilibrium  in  CC  A,  constant 
chemical  potentials  in  CC  B,  and  reasonable  transport  rates 
at  the  boundaries  between  the  CCs  and  the  membrane  region. 
Typically,  there  are  600  fluid  particles  in  CC  A  and  10-500 
fluid  particles  in  CC  B.  Simulations  were  started  with  par¬ 
ticles  distributed  in  CC  A  and  CC  B  only  and  no  particles  in 
the  membrane  region.  We  then  carried  out  typically  0.25  ns 
simulation  runs  to  achieve  steady  state.  Subsequent  produc¬ 
tion  runs  ranged  from  1  ns  for  pure  DCV-GCMD  simulations 
to  2  ns  for  mixture  DCV-GCMD  and  DCC-RxMD  simula¬ 
tions,  where  Af  =  3.438  fs. 

During  the  simulation  we  evaluated  the  excess  internal 
energy  u,  number  density  p,  pressure  P,  and  compositions 
xy  and  xf  of  the  mixtures  in  the  CCs.  The  pressure  was 
computed  from  the  virial  theorem.18  To  minimize  the  influ¬ 
ence  of  the  interfaces  between  the  CCs  and  the  membrane 
region  as  well  as  the  effects  of  the  confining  walls  on  the 
adjacent  portions  of  the  CCs,  fluid  properties  in  the  CCs 
were  calculated  in  a  predetermined  interior  portion  of  the 
CCs.  More  specifically,  the  fluid  property  calculations  were 
applied  only  to  molecules  in  CCs  whose  x  coordinates  were 
at  a  distance  greater  than  max{rc  ;R.}/= ,  from  the  CCs  bound¬ 
aries.  We  also  calculated  the  molar  flux  through  the  mem¬ 
branes  for  component  Z,  //,  and  the  permeability  of  compo¬ 


nent  Z,  Kt .  J i  was  determined  from  the  expression 

n\tk-n^ 

j  —  _ : _ 

1  NaA  yZN  MDsteps  A  t 

and  K ,  was  defined  as 


(34) 


Ki= 


Jx 

Ji 


(35) 


In  Eqs.  (34)  and  (35),  (Vj:l  K  and  /VfTL  are  the  net  movement 
of  particles  for  component  Z  through  the  yz  plane  at  x  =  0 
from  left-to-right  and  from  right-to-left,  respectively,  NA  is 
Avogadro’s  number,  Ayz  =  LxLz  is  the  yz  area  of  the  simula¬ 
tion  box,  ZVMDsteps  is  the  total  number  of  isothermal  MD 
steps,  and  A Pi=xfPA  —  xBPB  is  the  partial  pressure  differ¬ 
ence  for  component  Z. 


IV.  RESULTS  AND  DISCUSSIONS 
A.  DCV-GCMD 

Before  presenting  results  for  the  increase  of  reaction 
conversion  by  separation  of  the  H2  product  from  the  reaction 
mixture  of  CH4  /C02  /H2  /CO  for  the  dry  reforming  reaction, 
we  present  results  from  DCV-GCMD  simulations  for  the 
pure  fluids  and  equimolar  mixtures.  The  intent  is  to  investi¬ 
gate  the  properties  of  the  membrane  models  PM1,  PM2,  and 
RM  independent  of  chemical  reaction  behavior.  We  per¬ 
formed  all  simulations  at  a  typical  temperature  for  the  dry 
reforming  reaction,  T=  1100  K.20 

1.  Pure  fluids 

In  the  DCV-GCMD  simulations  for  pure  CH4,  C02, 
Hi,  and  CO,  we  set  p,A/(RT)  =  —4.039  and  varied  /jlb/(RT) 
from  —4.308  to  —5.520;  R  is  the  universal  gas  constant. 
This  results  in  a  pressure  difference  A  P  =  PA  —  PB  ranging 
from  7  bars  to  37  bars,  where  Pa^A%  bars.  Figure  4  presents 
J i  as  a  function  of  A P  for  all  three  membrane  models.  Values 
of  K i  for  the  membrane  models  did  not  show  (within  statis¬ 
tical  uncertainties)  dependence  on  AP.  Averaged  K /  values 
are  listed  in  Table  II.  Figure  4  and  Table  II  show  that  only 
PM1  is,  strictly  speaking  semipermeable,  i.e.,  permeable  to 
Hi  and  completely  impermeable  to  CH4,  C02,  and  CO. 
Both  PM2  and  RM  exhibit  very  small  yet  undesirable  per¬ 
meabilities  for  CH4 ,  COi ,  and  CO  in  addition  to  substantial 
permeability  for  H2 .  However,  semipermeable  membranes  in 
industrial  membrane  reactors  also  exhibit  small  undesirable 
permeabilities  for  nonseparating  components.  In  the  cases  of 
PM2  and  RM,  typical  values  of  the  molar  fluxes  for  CH4, 
C02 ,  and  CO  are  one  order  smaller  than  the  corresponding 
values  of  the  molar  flux  for  H2 .  Similar  behavior  is  exhibited 
for  component  permeabilities.  Note  that  since  crCH^crC02 
^ctc o  and  since  at  T=1100  K  values  of  e’s  have  only  a 
moderate  influence  on  the  transport  properties,  /ch4>-/co 
>Jqo2  in  the  case  of  PM2  is  primarily  due  to  mCH^<mC0 
<  in co,  ■  In  contrast  to  PM2  the  effect  of  the  different  masses 
of  the  components  on  J;  is  less  pronounced  for  RM.  This  is 
due  to  the  ability  of  the  particles  to  pass  more  freely  through 
the  straight  pores  of  PM2  as  opposed  to  the  contour  pores  of 
RM.  Also  note  in  Fig.  4  that  at  the  same  A P,  values  of  /Hi 
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FIG.  4.  Molar  flux  through  the  membrane  models  (a)  PM1,  (b)  PM2,  and  (c) 
RM  for  component  /,  7, ,  as  a  function  of  the  pressure  difference  A P  ob¬ 
tained  from  the  DCV-GCMD  simulations  for  pure  CH4  (■),  C02  (A),  H2 
(•),  and  CO  (♦).  Lines  serve  as  a  guide  to  the  eye  only. 


in  PM1  are  about  three  times  smaller  than  values  of  7Hi  in 
PM2.  Further,  values  of  /Hi  in  PM2  are  approximately  two 
times  larger  than  those  found  in  RM.  This  behavior  is  a  con¬ 
sequence  of  the  well  defined  pore  structure  of  PM2  whose 
PSD  is  much  narrower  compared  to  the  PSD  of  RM. 

Figure  5  shows  typical  density  profiles  for  H2  across  the 
membrane  models  PM1,  PM2,  and  RM  at  AP«=32  bars.  The 
high-pressure  and  low-pressure  sides  of  the  membrane  cor¬ 
respond  to  =®x  =  —  1  and  1  nm,  respectively,  implying  a  flux 
of  Ho  particles  from  left  to  right  in  Fig.  5.  The  vertical  dotted 
lines  denote  the  membrane  model  boundaries.  The  density 
profiles  for  PM1  and  PM2  exhibit  oscillatory  character  with 
minimums  corresponding  to  the  membrane  layers  and  with 
maximums  at  positions  between  them.  The  density  profile  for 
RM  decreases  nearly  monotonically  from  the  high-pressure 
to  low-pressure  membrane  sides  due  to  the  random  configu¬ 
ration  of  membrane  particles  in  RM.  The  behavior  of  the 
relative  permeabilities  of  H2  for  the  different  membrane 
models  is  reiterated  in  Fig.  5. 


Lisal  et  at. 

2.  Equimolar  mixtures 

In  the  DCV-GCMD  mixtures  simulations,  we  set 
{/if/(RT)}f=l={- 5.385,- 5.385,- 5.385,- 5.385}  and  var¬ 
ied  {/rf/(RT)}f=1  from  {-5.655,-5.655,-5.655,-5.655} 
to  {-6.732,-6.732,-6.732,-6.732}  by  keeping 
=  /J-(  ch  =  /if^  =  ju,Q0 .  This  results  in  AP  ranging  from  7  bars 
to  37  bars,  where  PA 50  bars.  At  T  =  1100  K,  mixtures  in 
CC  A  and  CC  B  are  marginally  nonideal  and  hence,  the  used 
values  of  /if  and  /if  produce  roughly  equimolar  mixture 
compositions  in  both  CCs,  i.e.,  xf^xf  =“0.25.  Plots  of  /,  as 
a  function  of  the  partial  pressure  difference  A Pt  for  all  three 
membrane  models  are  given  in  Fig.  6.  As  in  the  case  of  pure 
fluids,  values  of  K ,  did  not  show  (within  statistical  uncertain¬ 
ties)  dependence  on  A Pt  (see  Table  II).  From  Fig.  6  and 
Table  II  we  can  draw  similar  conclusions  about  the  depen¬ 
dence  of  the  molar  fluxes  and  permeabilities  based  on  the 
relative  membrane  PSDs  and  structures.  The  notable  differ¬ 
ence  being  that  the  values  of  7/  for  equimolar  mixtures  are 
smaller  than  corresponding  values  of  J t  for  pure  fluids  since 
AP;~0.25AP. 

B.  DCC-RxMD 

At  7’=1100  K  and  a  bulk  volume  V=VA  =  (L3x—L2x) 
XLVXL,=  1576.71  nm3,  both  the  REMC  and  RxMD  pre¬ 
dict  the  following  bulk  equilibrium  properties  of  the  reaction 
mixture  CH4 /C02 /H2 /CO:  u=  -0.05726  kJmol1,  P 
=  50.7n  bars,  p  =  0.3291184  nm~3,  xCh4  =  0.22210,  xCq2 
=  0.22210,  xy  =0.2781(),  and  a" ("o — 0. 278  ]  o  j  subscripts  in 
numerical  values  denote  the  standard  deviations  in  the  last 
digits.  Simulations  were  initiated  with  N'qq^= N'fff 

=N co=125  molecules  in  the  simulation  box.  At  equilib¬ 
rium,  the  total  number  of  molecules  in  the  bulk  reaction  sys¬ 
tem  was  /V=  ,/V/  =  5205. 

All  DCC-RxMD  simulations  were  started  with 
=Ni&=N%=N&=125  in  CC  A,  i.e.,  with  NT  =  500. 
Values  of  /if  were  chosen  in  such  a  way  to  obtain  virtually 
pure  H2  in  CC  B.  This  was  achieved  by  setting  /ifHJ(  KT) 
=  /ibcoJ{RT)  =  /ibcoI(RT)  =  - 12.111.  Values  of  /4iJ{RT) 
were  then  varied  from  —4.712  to  —8.078.  This  results  in 
AP>0  but  APH2<0  for  /ifj(RT)>- 5.25  and  AFy>0 
otherwise.  Note  that  the  pressure  in  the  reaction  void  de¬ 
creases  with  the  removal  of  H2  since  the  total  number  of 
particles  decreases  (see  Figs.  7-9).  For  APH2<0,  H2  flows 
from  the  transport  void  to  the  reaction  void  while  for  APHo 
>  0,  there  is  a  H2  flux  from  the  reaction  void  to  the  transport 
void.  Further,  /if/iRT)  =—  8.078  produces  nearly  vacuum 
conditions  in  CC  B,  i.e.,  the  equilibrium  number  of  particles 
in  CC  B  becomes  very  small  and  thus  PB—>0.  Hence,  the 
value  of  APHi  resulting  from  /if/iRT)  =  —  8.078  corre¬ 
sponds  roughly  to  the  maximum  achievable  APHi  . 

Figures  7-9  show  /Hi  and  NA/NAU ,  and  xf  as  a  func¬ 
tion  of  APHi  for  all  three  membrane  models.  Note  that  xf 
at  A  P h2  =  0  corresponds  to  xf  of  the  bulk  reaction  system. 
We  see  from  Figs.  7-9  that  the  maximal  APHi  [correspond¬ 
ing  to  /ift  /iRT)  =—  8.078]  is  directly  related  to  values  of 
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TABLE  II.  Permeabilities  K ,  obtained  from  the  DCV-GCMD  simulations  of  pure  CH4 ,  C02 ,  H2 ,  and  CO,  and  an  equimolar  mixture  of  CH4  /C02  / H2  /CO 
for  all  three  membranes  PM1,  PM2,  and  RM.  Subscripts  in  table  values  denote  the  standard  deviations  in  the  last  digits. 


Membrane 

Pure  fluids 

ZfCH4X106 

(mmol  bar-1  mm-1  s-1) 

/fco2X106 
(mmol  bar- 1  mm- 

‘s'1) 

/Th2X106 

(mmol  bar-1  mm-1  s-1) 

KcoX  106 

(mmol  bar  1  mm  1  s 

PMl 

09 

0, 

40050 

15u 

PM2 

20040 

7755 

127040 

145, 

RM 

35, 

3520 

57060 

35, 

Equimolar  mixture 

*ch4X106 

Kc  o2X106 

/tHixio6 

KcoX  106 

Membrane 

(mmol  bar-1  mm-1  s-1) 

(mmol  bar- 1  mm- 

‘s'1) 

(mmol  bar-1  mm-1  s-1) 

(mmol bar  'mm  1  s  ') 

PMl 

99 

0, 

45060 

0, 

PM2 

180, 

9525 

127070 

12040 

RM 

2530 

202o 

575 120 

3525 

/('ll  :  a  membrane  model  with  a  larger  /C[[i  has  a  lower 
A  1 1 1  (cf.,  e.g.,  the  case  of  PM2  with  1270X  10  6 

mmol  bar  1  mm-1  s  1  and  maximal  APHi— 3  bars  with  the 
case  of  PM1  with  K^— 400X  1CT6  mmolbar-1  mm'1  s_I 
and  maximal  APHo  —  8  bars).  Also  note  that  the  upper  por¬ 
tions  of  Figs.  7-9  show  that  the  values  of  /Hi  corresponding 
to  the  maximal  A PHi  are  approximately  the  same.  In  contrast 
to  these  values  of  /Ho ,  the  values  of  NA  IN'™  at  the  maximal 
A PHi  differ  significantly;  they  are  lower  for  membranes  with 
a  higher  K  u  .  Next  note  that  compositions  xj  as  a  function 
of  APHo  in  the  lower  portions  of  Figs.  7-9  show  that  the 
compositions  of  reactants  CH4  and  CCF  decrease  slowly 
with  increasing  APHo  (except  and  .tp0i  close  to  the 
maximal  A for  PM2).  The  composition  of  the  product 


FIG.  5.  Density  profiles  for  H2  across  the  membrane  models  PM1  ( - ), 

PM2  ( - ),  and  RNM  ( - )  at  APs«32  bars.  The  dotted  lines 

denote  the  membrane  boundaries.  H2  molecules  are  flowing  from  left-to- 
right  due  to  the  imposed  pressure  gradient. 


AP,  (bar) 


FIG.  6.  Molar  flux  through  the  membrane  models  (a)  PM1,  (b)  PM2,  and  (c) 
RM  for  component  /,  7/ ,  as  a  function  of  the  partial  pressure  difference  A P{ 
obtained  from  the  DCV-GCMD  simulations  for  equimolar  mixtures  of 
CH4(B)/C02(A)/H2(#)/C0(  ♦  ).  Lines  serve  as  a  guide  to  the  eye  only. 
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APh,  (bar) 


APh,  (bar) 


FIG.  7.  (a)  Hydrogen  molar  flux  ,/M  (•)  and  NA  INA'  (O)  as  a  function  of 
the  hydrogen  partial  pressure  difference  APHi,  and  (b)  the  composition  in 
the  reaction  void  xf  (CH4 ,  ■;  H2 ,  •;  CO,  ♦ )  as  a  function  of  APH  in  the 
case  of  PM1  obtained  from  the  DCC-RxMD  simulations.  Due  to  identical 
initial  compositions  of  CH4  and  C02  x cn4 = *co2  within  statistical  uncertain¬ 
ties,  therefore  Xq0  is  not  plotted.  Lines  serve  as  a  guide  to  the  eye  only. 

CO  increases  significantly  with  increasing  A.PH2 .  At  the 
maximal  A.PHi  ,  Xqq  is  larger  for  membranes  with  a  larger 
/('ll  .  With  respect  to  x'^0  at  APHo  =  0,  increases  of  x^q  at  the 
maximal  A  PH^  are  —60%  for  PM1,  —80%  for  RM,  and 
— 150%  for  PM2.  x^4  decreases  with  increasing  APHi  due  to 
hydrogen  separation  from  the  reaction  void  to  the  transport 
void.  The  decrease  of  is  related  to  K  and  is  larger  for 
membrane  models  with  larger  Ku  .  The  total  yield  of  H2  is  a 
result  of  the  H2  amount  in  both  the  reaction  and  transport 
voids. 

V.  CONCLUSIONS 

We  have  presented  a  simulation  tool,  termed  the  dual 
control  cell  reaction  ensemble  molecular  dynamics  (DCC- 
RxMD)  method,  to  study  fluid  mixtures  that  are  simulta¬ 
neously  chemically  reacting  and  adsorbing  in  a  porous  ma¬ 
terial.  The  DCC-RxMD  method  was  developed  by  coupling 
a  nonequilibrium  molecular  dynamics  method  with  two 
Monte  Carlo  based  methods,  namely,  reaction  ensemble 
Monte  Carlo  (REMC)  and  grand  canonical  Monte  Carlo 
(GCMC).  Control  cells,  which  are  in  direct  physical  contact 
with  the  porous  solid,  are  used  to  maintain  the  desired  reac- 


Lisal  et  at. 


APH;  (bar) 


aPHj  (bar) 


FIG.  8.  (a)  Hydrogen  molar  flux  Ju  (•)  and  N A  INA"  (O)  as  a  function  of 
the  hydrogen  partial  pressure  difference  APHi,  and  (b)  the  composition  in 
the  reaction  void  xf  (CH4 ,  H2 ,  CO,  ♦  )  as  a  function  of  APH  in  the 
case  of  PM2  obtained  from  the  DCC-RxMD  simulations.  Due  to  identical 
initial  compositions  of  CH4  and  C02  a('i2  =  xf0^  within  statistical  uncertain¬ 
ties,  therefore  xf()  is  not  plotted.  Lines  serve  as  a  guide  to  the  eye  only. 

tion  and  flow  conditions.  The  simulation  setup  closely  mim¬ 
ics  an  actual  experimental  system  in  which  the  thermody¬ 
namic  and  flow  parameters  are  precisely  controlled.  The 
method  is  akin  to  the  dual  control  volume  grand  canonical 
molecular  dynamics  method  that  was  developed  to  study  the 
transport  properties  of  fluid  mixtures  primarily  in  confined 
systems.  The  added  feature  of  the  DCC-RxMD  method  is  the 
inclusion  of  chemical  reactions,  thus  its  applicability  is  to  a 
wider  range  of  processes  beyond  solely  adsorption  phenom¬ 
ena.  The  method  presented  here  allows  for  the  calculation  of 
both  equilibrium  and  nonequilibrium  transport  properties  in 
porous  materials  such  as  diffusion  coefficients,  permeability, 
and  mass  flux.  Effects  on  these  properties  due  to  the  charac¬ 
teristics  of  the  porous  material  can  be  predicted;  characteris¬ 
tics  such  as  the  pore  size  distribution,  connectivity,  porosity, 
and  surface  area.  Note  that  the  DCC-RxMD  method  reduces 
to  the  RxMD  method11  if  the  REMC  steps  are  performed  in 
both  control  cells.  Also  note  that  in  general,  neither  the 
DCC-RxMD  nor  the  RxMD  methods  can  provide  reaction 
rate  information  but  in  turn,  neither  method  is  limited  by 
reaction  rates  or  activation  energy  barriers. 

As  an  illustration  of  the  method,  we  simulated  the  dry 
reforming  of  methane  reaction  within  a  nanoscale  reactor 
model  in  the  presence  of  a  semipermeable  membrane  that 
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AP„;  (bar) 


aPHj  (bar) 


FIG.  9.  (a)  Hydrogen  molar  flux  ,/M  (•)  and  NA  IN'f  (O)  as  a  function  of 
the  hydrogen  partial  pressure  difference  A  /’  Hi ,  and  (b)  the  composition  in 
the  reaction  void  xf  (CH4 ,  ■;  H2 ,  •;  CO,  ♦ )  as  a  function  of  A  Pu  in  the 
case  of  RM  obtained  from  the  DCC-RxMD  simulations.  Due  to  identical 
initial  compositions  of  CH4  and  C02  Xq^=x/(-'.()_  within  statistical  uncertain¬ 
ties,  therefore  Xq0  is  not  plotted.  Lines  serve  as  a  guide  to  the  eye  only. 

was  modeled  as  silicalite.  We  studied  the  effects  of  the  mem¬ 
brane  structure  and  porosity  on  the  reaction  species  perme¬ 
ability  by  considering  three  different  membrane  models.  We 
also  studied  the  effects  of  an  imposed  pressure  gradient 
across  the  membrane  on  the  mass  flux  of  the  reaction  spe¬ 
cies.  The  conversion  of  syngas  (H2/CO)  increased  signifi¬ 
cantly  in  all  the  nanomembrane  reactor  models  considered. 
An  increase  of  60%  at  the  maximal  imposed  pressure  oc¬ 
curred  in  the  nanomembrane  reactor  model  with  a  membrane 
that  was  exclusively  permeable  to  H2 .  The  larger  increase 
(80%  and  150%  at  the  maximal  imposed  pressures)  occurred 
in  the  nanomembrane  reactor  models  with  membranes  that 
exhibited  larger  values  of  H2  permeability  but  at  the  same 
time  possessed  very  small  yet  undesirable  permeabilities  for 
the  nonseparating  components  CH4 ,  C02 ,  and  CO. 

In  the  application  of  the  DCC-RxMD  method  presented 
here,  we  considered  membrane  models  that  did  not  exhibit 
attractive  character  towards  the  reaction  species.  Further 
simulations  that  include  effects  of  the  relative  adsorption  be¬ 
havior  of  the  reaction  species  would  be  worthwhile.  Here 
again  such  simulations  could  provide  insight  into  the  mem¬ 
brane  characteristics  that  would  most  influence  the  conver¬ 
sion  of  reactants  to  products.  Further,  we  note  that  analogous 
to  the  REMC  method,  multiple  reactions  can  be  simulated 


simultaneously  in  the  DCC-RxMD  method;  providing  a 
means  of  studying  systems  in  which  competing  reactions 
play  a  crucial  role.  Finally,  we  should  note  that  in  this  study, 
the  DCC-RxMD  method  was  applied  to  a  system  where  a 
semipermeable  membrane  exists  between  the  reaction  and 
transport  voids.  However,  this  could  be  any  type  of  porous 
solid  model,  be  it  a  simpler  slitpore  model21  or  a  more  real¬ 
istic  model  such  as  an  activated  carbon.40  Within  the  limit  of 
the  available  computational  resources,  DCC-RxMD  simula¬ 
tions  could  also  include  porous  solids  on  a  longer  length 
scale  such  as  microporous  solids. 
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